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Abstract 

Random boolean networks are a model of genetic regulatory networks that has proven 
able to describe experimental data in biology. They not only reproduce important phe¬ 
nomena in cell dynamics, but they are also extremely interesting from a theoretical view¬ 
point, since it is possible to tune their asymptotic behaviour from order to disorder. The 
usual approach characterizes network families as a whole, either by means of static or 
dynamic measures. We show here that a more detailed study, based on the properties of 
system’s attractors, can provide information that makes it possible to predict with higher 
precision important properties, such as system’s response to gene knock-out. A new set 
of principled measures is introduced, that explains some puzzling behaviours of these net¬ 
works. These results are not limited to random Boolean network models, but they are 
general and hold for any discrete model exhibiting similar dynamical characteristics. 


1 Introduction 

Biological systems typically involve a great number of interacting units exhibiting a high de¬ 
gree of self-organisation that ensures their continuous functioning and allows them to adap¬ 
tively respond to environmental changes. Understanding the roots of such nontrivial proper¬ 
ties represents a major challenge in theoretical biology and complex systems science. Among 
these systems, a prominent role is that of genetic regulatory networks, which in multicellular 
beings rule the expression (i.e., the transcription) of the various genes. A major question 
related to GRNs behaviour is the coexistence, in multicellular organisms, of several different 
cell types with a single common genome. It has been proposed [mEO] that this could be 
explained by the existence of several possible dynamical behaviours in the same dynamical 
system. In particular, it has been suggested that different cell types correspond to different 
asymptotic dynamical patterns of activity (attractors) of the same system. According to this 
view, the genome dictates the units and their basic mutual interactions, and the attractors de¬ 
termine the cell type. Note that this picture can be expanded to accommodate also epigenetic 
effects. 
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A prominent model of genetic regulation is that of Random Boolean networks (RBNs), 
which were introduced by Kauffman imiiB and became one of the major models of complex 
systems due to their interesting dynamical behaviour [HI EQl El H El E]. The interest later 
faded, but in recent years it has been renewed by important theoretical advances [MlEZlIlslIIll 
and also, as far as the application to genetics is concerned, by the availability of genome-wide 
expression data which can be properly described by RBNs [Ml[MlEH EH ESIIT]. Further 
reasons of this renewed interest are related to the possibility to describe complex phenomena, 
such as cell differentiation |42l [28l 00] , and other biological systems such as whole organisms 
or tissues [MiEHiEniiinii. 

RBNs support different dynamical regimes and are able to describe the transition from 
ordered to disordered phases by changing the values of some of their parameters. R is well 
known that the bias (i.e., the internal homogeneity) and the average number of input gene 
variables (i.e., the network connectivity) can modulate the order/disorder transition, with 
higher homogeneity and lower connectivity leading toward more ordered behaviour. In ad¬ 
dition, also the choice of the functions associated to each node plays a very important role; 
indeed, RBN dynamical behaviours are deeply influenced by the various classes of the Boolean 
functions chosen (see Section for more details). It has become commonplace to characterize 
RBNs as either ordered, critical or disordered depending upon their asymptotic dynamics. 
Typically, the attractors of ordered networks are stable, in the sense that they can be reached 
by nearby initial conditions, while disordered networks display a behaviour that is the dis¬ 
crete analogue of the butterfly effect, i.e., a small perturbation of the initial condition usually 
leads to a different attractor. Because of this latter property, disordered networks are also 
called chaotic, with slight abuse of term. Critical networks display intermediate behaviours 
and are often considered the ones which can best capture important properties in biological 
systems. However, due to their large degree of randomness, the dynamical characterization 
of these systems is a subtle issue. Commonly, they are characterized on the basis of their 
structural parameters, thereby providing what is actually a characterization of a whole family 
(or ensemble) of such networks. But single network realizations may behave in a way very 
different from the one that is typical of the family. Moreover, even single attractors of the 
same network can behave in a way that is different from the others. Obviously, the dynamical 
behaviours of the attractors depend on system’s structure, but a common origin does not 
imply that these dynamical behaviours have to be to similar. 

In this paper, we introduce a detailed way to describe the behaviour of RBNs and we will 
clarify the relationship between our results and those obtained by the usual ensemble analy¬ 
sis. We will show that this new set of measures is of great help in clarifying and explaining 
situations that cannot be properly dealt with the usual methods. 

The paper is organized as follows: Section [^briefly presents the RBN framework; Sectionj^ 
presents the more common static and dynamic measures of the dynamical RBN regimes, which 
are deeply discussed in Section [^ Section presents our detailed vision of the behaviours 
observed in the so-called critical RBNs. In Section we discuss a particularly significant case 
in which the new measures allow a new interpretation of experimental biological tests, and 
Section [3 summarises and concludes the work. 
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2 Random Boolean network models 


In this section we succinctly review RBNs, emphasizing the most relevant properties and 
results for this work. The interested reader is referred to existing reviews for a more thorough 
discussion of RBNs [iziiiiiisniiaEi. 

A Boolean network (BN) is a discrete-state and discrete-time dynamical system whose 
structure is defined by a directed graph of N nodes, each associated to a Boolean variable Xi, 
i = 1,..., N, and a Boolean function FAxi ,,..., x*. ), where kin i is the number of inputs of 

node i. The arguments of the Boolean function Fi are the values of the nodes whose outgoing 
arcs are connected to node i. The state of the system at time t, f G N, is defined by the array 
of the N Boolean variable values at time t. s{t) = (xi(t),... ,XAr(f)). The most studied BN 
models are characterized by having a synchronous dynamics—i.e., nodes update their states 
at the same instant—and deterministic functions. However, many variants exist, including 
asynchronous and probabilistic update rules [35] . 

A special category of BNs that has received particular attention is that of RBNs, which 
can capture relevant phenomena in genetic and cellular mechanisms and complex systems in 
general. RBNs are usually generated by choosing at random kin inputs per node and by defin¬ 
ing the Boolean functions by assigning to each entry of the truth tables a 1 with probability p 
and a 0 with probability 1 — p. The parameter p is referred to as the bias. Depending on the 
values of kin and p the dynamics of RBNs is called either ordered or ehaotic. In the first case, 
the majority of nodes in the attractor is frozen and any moderate-size perturbation is rapidly 
dampened and the network returns to its original attractor. Conversely, in chaotic dynamics, 
attractor cycles are very long and the system is extremely sensitive to small perturbations: 
slightly different initial states lead to divergent trajectories in the state space. RBNs temporal 
evolution undergo a second order phase transition between order and chaos, governed by the 
following relation between kin and p: 

= [2pc(l(1) 
where the superscript c denotes the critical values m- 

3 The measure of RBN dynamical regimes 

In this section we hrst introduce and discuss the usual ways for characterizing the dynamical 
regimes of RBNs. Then, we revise the notion of critical network showing that the dynamic 
behaviour of a BN can be dramatically different across its attractors. 

3.1 The spread of perturbations 

As previously observed, ordered and disordered dynamical regimes are usually described as 
systems’ behaviours supporting respectively short attractors with fairly regular basins of 
attraetion^ and long attractors with sensitive dependence upon initial conditions 0 There 
are two different ways to identify ordered and disordered regions: (i) a static one, based upon 

^An attractor basin is the set of states whose evolution lead to the attractor, its size (or dimension) being 
the cardinality of the set. 

^We remind that in disordered systems trajectories starting from nearby points typically lead to different 
attractors. 
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the knowledge of the nodes average connectivity and the bias of the Boolean functions and 
(a) a dynamical one based upon the study of the spreading of perturbations through the 
system. 

In ordered networks a change at one node (i.e., a bit flip) propagates in one step on 
average to less than one other node: starting from a random initial condition an ordered 
system rapidly reaches a stable condition in which the majority of nodes are frozen; if this 
asymptotic behaviour is perturbed, there is a very high probability of coming back to it. On 
the contrary, in disordered networks a perturbation at one node propagates in one time step 
on average to more than one other node: very close initial conditions could rapidly diverge 
toward different attractors, and attractors typically have long periods, with large portions 
of oscillating nodes; if perturbed, the system has therefore high probability of changing its 
asymptotic behaviour. Critical systems are at the boundary between these two dynamical 
regimes: a change at one node propagates in one time step on an average to exactly one other 
node. This situation corresponds to the percolation of frozen nodes through the network, and 
therefore to the formation of wide but isolated oscillating zones. 

3.2 Static vs. dynamic estimates 

The main static methods to measure the RBN dynamical regimes implicitly presume ergod- 
icity, that is, inputs arise with the same probability during evolution, and time average over 
the states visited by the network yields the same results as averages over the whole phase 
space. A widely used measure is the average sensitivity of a Boolean function, proposed by 
Shmulevich and Kauffman m- The sensitivity of a function T) measures how sensitive the 
output of the function is to changes of its inputs. Let us consider an input configuration x and 
all its 1-Hamming neighbours, i.e., the input configurations that differ from x in exactly one 
value; the sensitivity s^'{x) is defined as the number of 1-Hamming neighbours of x on which 
the Fi function values are different than on x 137|. The average sensitivity of the function T), 
I{Fi) is the expected value of sensitivity with respect to the distribution of x. Related 

to the average sensitivity is the notion of influence of a variable: the influence of the j-th 
input variable of a function T), Ij{Fi), is the probability that the function T) changes its value 
when the value of the j-th variable is changed (a concept linked to Boolean derivatives and 
Lyapunov exponents of RBNs 0 El]). It can be proven that /(T)) is kin times the average 
probability that the output of a node changes when one of its inputs changes, in formulas: 

I{Fi) = (1 - qi)kin,i (2) 

where qi is the probability that the output of the i — th node does not change when one of 
its kin^i inputs changes m- 

At this point Shmulevich and Kauffman introduce the ergodic hypothesis: under uniform 
input distribution the average sensitivity of function T) is simply equal to the sum of the 
influences of all its input variables (a number that ranges from 0 to kin)- The average 
sensitivity of a BN is the weighted average of the sensitivities of all its functions. If this final 
average sensitivity is lower, equal to or higher than 1 the network typically dampens, maintains 
or amplifies the perturbations and is therefore ordered, critical or disordered, respectively. 
It should be noted that, when the Boolean functions are drawn according to a Bernoulli 
distribution with a specific bias p, the average sensitivity coincides with the well-known critical 
transition curve kin{‘2.p{l — p)) (see |37j). However, when network functions are generated 
according to probability distributions that favour some variables relative to others, or when 
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functions are chosen randomly from certain classes of functions (e.g., canalizing), the average 
sensitivity well captures the dynamical behaviour of the network. We also remind that the 
ergodicity assumption does not hold for the dynamics of arbitrary RBNs |23[ I39j . 

The alternative to static measures is that of explicitly exploring the dynamical behaviour of 
the system. An interesting and well-known method is the so-called Derrida procedure mm 
which exploits the spreading of perturbations through the network. This procedure involves 
two parallel rnns of the same system, with initial states that differ for only a small fraction 
of the nnits. This difference is usually measured by means of the Hamming distance. Let 
h{t) be the Hamming distance between the states at time t in the two parallel trajectories. 
If after a transient the two runs are likely to converge to the same statej^ i.e., h{t) 0, 
then the dynamics of the system is robust with respect to small perturbations—a signature 
of the ordered regime. Conversely, if the difference is likely to increase, then the dynamics 
is sensitive to small perturbations and the corresponding regime is disordered. The common 
way to measure the dynamical regime of a RBN by means of the Derrida procedure is that 
of randomly generating a great number of pairs of initial conditions differing for one or a 
few units, evolving the network for one step, measnring the Hamming distance of the two 
resulting states, taking the averages for each perturbation size, and plotting these averages 
(i.e., h{l)) vs. the initial perturbation size (i.e., h{0)). Let A be the slope at the origin of the 
cnrve /i(l) vs. h{0). A system is ordered for A < 1, whereas it is disordered for A > 1. The 
case of A = 1 identifies the critical regime [2] . It is worth stressing that the Derrida parameter 
is an empirical estimation of the average sensitivity mentioned above. Furthermore, under 
the hypothesis of ergodicity, it holds A ~ (1 — q)kin, where q is the average probability that 
the output of a node does not change when one of its inputs does m- 

3.3 Attractor sensitivity 

It is important to notice that the usual interpretation of the behaviour of RBNs, e.g., as 
models of genetic regulatory networks, privileges their attractor states. Indeed, the network 
will be found in one of these states, after transients have died out. Therefore, measures 
taken on randomly chosen states do not necessarily provide a correct estimate of the relevant 
dynamical behaviour of the system, which should rather be evaluated on its attractors. 

A more meaningful way to analyse the dynamics is that of applying the classical Derrida 
procedure only on the states belonging to the attractors um- We can therefore define the 
sensitivity on attractor i (S'A) as the result of the Derrida procedure performed only on the 
states belonging to the attractor i, and the sensitivity on attractors {SA) as the average of 
the SAi, each SAi being weighted with the size of its attraction basin. We will also refer to 
the usual Derrida procedure, i.e., the one on random initial states, as DA. 

A similar approach was used in |15j . where the anthors only considered what has been 
defined here as the quantity SA. As we will see in the following, the same RBN can support 
several attractors with different dynamical behaviour, therefore the SA value alone does not 
provide a complete picture and, in some conditions, can even snggest misleading conclusions 
on the dynamical regime. In order to appreciate differences and similarities among DA, 
SA and SAi, we will analyse these measnres across different statistical ensembles. In the 
following we use the word family for networks with transition functions randomly picked from 
the same set J-. Families will be denoted by the symbols AIi, AI 2 ) • ■ ■ j Adi- We assume that 

^The estimation is performed on many different initial conditions. 


5 



1.4 



Classic Derrida 


Figure 1: Attractor sensitivity of the attractors of 40 networks. The networks are grouped in 
4 sets having different number of nodes. As it can be noted, the attractor sensitivities of nets 
having the same DA (vertical sequences of marks) can have very different values. 


Table 1: The averages of DA and of SA for the 4 sets of nets in Figure and 


N 

{DA) 

(5A) 

50 

1.00T0.09 

l.OTO.l 

100 

0.99T0.06 

1.00T0.05 

150 

1.00T0.05 

1.01T0.05 

200 

1.00T0.05 

1.01T0.05 


each node has the same number of incoming links kin, and that the origin of these links is 
chosen at random with uniform probability among the remaining N — 1 nodes, prohibiting 
multiple connections. As previously mentioned, in random networks by varying D and/or p 
it is possible to move from disordered to ordered dynamical regimes. 

Let us first consider the usual ensemble, i.e., families composed of networks of different 
size with random topology, kin = 2 and Boolean functions uniformly distribnted {p = 0.5), 
referred to as A4i family in the following. As one can observe in Figure[^ the attractors of the 
same network can have values of SAi that significantly differ from one another and that in turn 
may be very different from the corresponding DA value of the network. These differences can 
be observed in several networks across different families characterized by different number 
of nodes. However, one can also see from Figure that SA (i.e. the weighted average 
of the attractors’ sensitivities) is strongly correlated with DA (note that the relation is well 
approximated by a linear function of slope 1). This result means that the dynamical behaviour 
averaged across all the attractors of a single RBN is related to the dynamical behaviour of a 
random sampling of the systems’ state space. 

Even more remarkably, from the results in Table we observe that by averaging over 
networks with the same size, the average values of DA and of S'A are almost equal|^ A 

^In this case, the averages are around 1, as the ensemble considered is the first historical example of critical 
systems. 
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Figure 2: The SA values of the same nets of Figure SA correlates with DA so that there 
is a clear—although noisy—proportionality between SA and DA. The figure shows also the 
bisector (the continuous line). 


similar situation holds for other families generated with different values of kin and p, in all 
the three dynamical regimes. In all these families, the attractor sensitivity considerably differs 
even across attractors of the same network. However, DA ~ SA when the values are averaged 
across all the networks of the family. 

It is important to observe that this last property does not hold for all the possible fam¬ 
ilies. Let us now consider two families of RBNs—in the following denoted by AI 5 and Me, 
respectively—with kin = 2, in which only a subset of canalizing Boolean functions is allowed. 
The interest for these particular families comes from a study on random threshold networks 
(RTNs), in which each node computes the weighted sum of its input (we consider the case of 
weights in { — 1 , +!})• If the sum exceeds a given threshold, the node takes the value 1 ; 0 , oth¬ 
erwise. It is possible to map these transition rules on the RBN framework, obtaining a set J- of 
rules for each threshold. In this work, we study the dynamical behaviour of the RBN families 
M 5 and Me , whose activation patterns correspond to RTNs having respectively thresholds 
h = -1-0.5 and h = —0.5. Incidentally, the two sets of Boolean functions identified in this way 
are complementary to each other. Despite this similarity, the asymptotic behaviours of these 
two families are quite different, and their SA value is considerably different from the classical 
Derrida measure—see Table for the details. Indeed, the attractor sensitivity is 0.96 for the 
Me family and 0.65 for the Me family (70 nodes per network). This difference can provide 
an explanation for the striking different dynamical behaviour in the two families: the average 
number of the attractors of networks in Me is significantly higher than the average number of 
attractors \n Me] considerable differences can also be found in length of cycles and number of 
frozen nodes [ 8 ]. The two families of RBNs show therefore very different behaviours despite 
the fact that both have the same DA. Hence, the Derrida parameter is not able to correctly 
describe the dynamics in these cases, where there is a remarkable difference between DA and 
SA. 

An even more extreme case in which DA and SA are very different concerns the case of BNs 
evolved to perform a global computation task. We consider the case of the so-called Density 
Classification Problem, which consists of classifying a binary string in either of two classes. 
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Table 2: The Boolean functions allowed in the and families {N = 70, kin = 2) and 
their respective measures of dynamical behaviours. 



M 5 

Me 

A B 

Ay B ^Aab Aa^b false 

-.(TVS) -.TVS AV-.B TRUE 

0 0 

0 0 0 0 

1111 

0 1 

110 0 

0 10 1 

1 0 

10 10 

0 0 11 

1 1 

10 0 0 

0 111 

average sensitivitj ^ 

0.75 

0.75 

DA 

0.74 

0.74 

SA 

0.96 

0.67 


depending on the ratio between Os and Is. The goal is then to find the Boolean functions of a 
network such that it is driven to a uniform state, consisting of all Is, if the initial configuration 
contains more Is and all Os otherwise [23]. It can be shown that the simple majority rule 
applied on random topologies outperforms all human or artificially-evolved rules running on 
ordered lattice, with a performance that increases with the size of the network [271122]. The 
majority rule states that the value of a node at time f -|- 1 is 0 (resp. 1) if the majority of its 
neighbours has value 0 (resp. 1) at time t. In this context, we studied a family of BNs (Mr), 
with kin = 3 and N = 71evolved to solve this task (see |6] for details). We are interested 
here in the analysis of the evolved networks in terms of sensitivity and Derrida parameter. 
According to the value of DA, this family is chaotic (with DA = 1.50), but the attractor 
sensitivity SA supports a very different conclusion: we have SA < 0.001, indicating that the 
system is deeply in the ordered region. And indeed the system is very ordered, having just a 
few very short attractors (typically, only two fixed points) with regular basins of attraction 
(nearby initial conditions typically evolve towards the same attractor). Figuredepicts the 
distributions of DA and SA in 200 realizations, showing the striking difference existing in 
this case between these parameters. 

In summary, we remark that in general DA and SA can considerably differ. Furthermore, 
even when they coincide, we must be very careful in drawing conclusions on the dynamical 
regimes, as the individual sensitivities on single attractors can be very dissimilar. In the next 
section we discuss the reasons for these differences and we outline their implications. 

4 Linking static and dynamical measures 

The point where static and dynamical measures deviate is that of the assumption of ergodicity 
made by the static approach. As previously pointed out, under the hypothesis of uniform 
input distribution, it can be shown that the average sensitivity of a function F) is equal to the 
sum of the influences of all its input variables ISZj. Nevertheless, in all the other cases this 
computation requires more attention. In particular, the asymptotic states sample only a part 
of the whole state space, a part that in non-ergodic systems can be quite limited. In order 
to circumvent this problem one needs a way to correctly weight each input conhguration. In 
RBNs it is possible to estimate the SAi of the f-th attractor from its time series by computing 
the average fraction of “1” (i.e., its frequency of occurrence h), and weight in such a way the 
influence of the input variables of each function. To obtain the average sensitivity of function 

®An odd number of nodes avoids the cases with equal quantity of Os and Is 
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Figure 3: Probability distribution of classical Derrida values (DA) and of attractor sensitivity 
(5^) for RBNs of 71 nodes in the family Aij. 200 networks sampled. 


Fi, i.e., I{Fi), we sum these weighted influence values and then we take the global average over 
all the Boolean functions, weighting each function by its occurrence in the network (Table [^. 
The global measure SA is computed by averaging over all the attractor sensitivity values, 
SAi, weighting each SAi with the basin of attraction size of the corresponding attractor. 
Table and Table show the values of a particular SAi for and Mq families, whereas 
Table shows the dependence of SA from the number of nodes N. 

We observe that b, i.e., the asymptotic number of “1”, depends on the dynamics of the 
system, which in turn originates from system’s features such as connectivity and Boolean 
functions. Therefore, it should be possible to analytically estimate this value. Indeed, one may 
use the so-called annealed approximation, which is a sort of mean field approximation—that 
holds for annealed networks HU with an infinite number of nodes—to estimate b and hence 
the value of SA. This approach can give reasonable guesses also for quenched systems [TTl 
I39| . Nevertheless, it takes averages and therefore ignores the possible peculiar behaviour on 
individual attractors; moreover, it cannot take into account the effect of the finite number 
of nodes that may be often non-negligible. An example of this discrepancy can be observed 
for the family (see Table [^. Eventually, by knowing the asymptotic sensitivity I{Fj) of 
function Fj and by using Equation it is possible to derive for each node j the probability of 
output change in case of an input change (or the complementary probability qj of maintaining 
its value constant). 

5 “Critical” networks can have very different dynamical be¬ 
haviours 

The aforementioned considerations provide us with instruments to appreciate in finer de¬ 
tails the differences among the networks usually classified “at the edge of chaos” [19]: the 
values of their static parameters are all “critical”, but the dynamical behaviours of their 
attractors could significantly deviate from critical regimes. This phenomenon finds striking 
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Table 3; The table shows the probability of occurrence P{A, B) of each input configuration 
vs. b, i.e., the normalized occurrence of value 1. A “CA” is reported if a change on the flipped 


input has as a consequence 

a change 

in the correspond 

[ing function output. 

Input 

P{A,B) 

Family At 5 

Family Aie 

A* 

B 


Ay B 

A A —iB 

“lA A B 

FALSE 

“'(A V B) 

—'A V B 

A V ^B 

TRUE 

0 

0 

{i-br 

CA 

CA 

nc 

nc 

CA 

CA 

nc 

nc 

0 

1 

6(1 - 6 ) 

nc 

nc 

CA 

nc 

nc 

nc 

CA 

nc 

1 

0 

6(1 - 6 ) 

CA 

CA 

nc 

nc 

CA 

CA 

nc 

nc 

1 

1 

6 ^ 

nc 

nc 

CA 

nc 

nc 

nc 

CA 

nc 

A* 

B 


Ay B 

A A ^B 

-^A A B 

FALSE 

“■(A V B) 

—'A V B 

A V ^B 

TRUE 

0 

0 

( 1 - 6 )^ 

CA 

nc 

CA 

nc 

CA 

CA 

CA 

nc 

0 

1 

6(1 - 6 ) 

CA 

nc 

CA 

nc 

CA 

nc 

CA 

nc 

1 

0 

6(1 - 6 ) 

nc 

CA 

nc 

nc 

nc 

CA 

nc 

nc 

1 

1 


nc 

CA 

nc 

nc 

nc 

CA 

nc 

nc 


Table 4: The table shows the influence of the input variables, computed by taking into 
account the effective occurrence probability of the possible input configurations, and the 
resultant function sensitivity. The sensitivity of the considered attractor, SAi, is estimated (by 
construction, in each family the four Boolean functions have the same occurrence probability). 



Family M 5 (6 = 0.08) 

Family Me (6 = 0.67) 


Ay B 

AA^B 

-A A 5 

FALSE 

AAy B) 

-.A VB 

Ay ^B 

TRUE 

Influence of A 

0.92 

0.08 

0.92 

0 

0.33 

0.22 

0.67 

0 

Influence of B 

0.92 

0.92 

0.08 

0 

0.33 

0.78 

0.33 

0 

Function sensitivity: 

1.84 

1 

1 

0 

0.66 

1 

1 

0 

SAf. 

0.96 




0.67 





Table 5: For each family the table shows the theoretical estimate of SA (computed by means 
of the procedure explained in the text) and the corresponding experimental value (estimated 
by effectively performing the Derrida procedure on random initial conditions and only on the 
attractor states respectively for DA and SA measures), vs. different number of nodes. 


Measure 

N 

b 

Family AA 5 

b 

Family A4q 




Theo. 

Exper. 


Theo. 

Exper. 

DA 

70 

0.5 

0.75 

0.74 

0.5 

0.75 

0.74 

SA 

70 

0.08 

0.96 

0.93 

0.67 

0.66 

0.64 

SA 

700 

0.05 

0.97 

0.95 

0.66 

0.67 

0.67 

SA 

00 

0 

1 

— 

0.67 

0.67 

— 
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(c) 


Figure 4: Measures on “critical” families (their bias p deriving from Equation]^ with N = 200 
and kin variable from 2 to 10, 50 networks for each family, (a) Number of attractors in 
each RBN (average, median, maximum); (b) median values of DA and SA] (c) minimum 
and maximum value of DA and SA] (d) the dispersion (standard deviation) of the same 
distribution in (b) and (c). 
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evidence in the results of a series of experiment we made, in which we created 50 networks 
for kin £ {2,3,..., 10} and N = 200, with Boolean functions generated by using the bias 
p = Pc, where pc is the critical value derived from Equation [Tj In Figure [^, we can see 
that the statistical ensembles show a critical behaviour, as the measures DA and SA assume 
values very close to 1. Nevertheless, the number of attractors in each RBN class dramatically 
decreases as kin increases (Figure]^). Furthermore, the minimum and maximum values of 
DA remain constant as kin increases, whereas the corresponding minimum and maximum 
values of SA increase their (already remarkable) distance from 1.0 (see Figure]^); this fact 
can also be observed looking at the dispersion of their respective distributions (the standard 
deviations of DA and distributions shown in Figure]^). Obviously, the single attractor 
sensitivities SAi are still more extreme . Therefore, although the statistical ensembles main¬ 
tain a critical behaviour across the values of kin, each single RBN realization has higher and 
higher probability to deviate from the critical behaviour; furthermore, the magnitude of this 
deviation grows as kin increases. Interestingly, the same deviation increase can be observed 
also in RBNs whose bias p differs from the “edge of chaos” values; for example, we observed 
ordered RBN ensembles in which, by increasing kin, one observes that there are progressively 
more and more RBN realizations with differing DA and SA, although their ensemble aver¬ 
ages remain close to each other. On the basis of these results, we can claim that by using 
these measures it is possible to introduce a finer classification within the “classical” critical 
networks that allows a clearer understanding of their asymptotic dynamical behaviour. In 
addition, it is worth to emphasize that the discrepancy between this finer analysis and the 
usual one increases as kin increases. 

6 Consequences for biological behaviours 

The previous considerations have important consequences on the application of RBNs as 
model of genetic regulatory networks. It is very interesting that, on the same network real¬ 
ization, DA and SA can have different values, but even more remarkable is the fact that the 
same RBN can support attractors showing different attractor sensitivities and therefore dif¬ 
ferent dynamical behaviours, that can hopefully be accessible to experimental measures. As 
a leading example, we can use a prominent interpretation of biological measures, performed 
by using the RBN framework to describe gene knock-out experiments in Saccharomyces cere- 
visiae [Ml ESI El]. 

A gene knock-out is a permanent silencing of a single gene: biologists can experimen¬ 
tally induce such modification and quantify (by using for example cDNA microarrays) the 
corresponding change of expression of all the other genes (the so-called avalanche) [16| . It is 
possible to reproduce in silica the knock-out process: at a certain time point, starting from 
a state of an attractorQ the value of one of the nodes of a RBN is permanently clamped to 
the value 0. The evolution of the unperturbed network is compared to that of the perturbed 
one (after the latter has reached again a stable asymptotic behaviour). A gene is said to be 
affected (or perturbed) if its behaviour differs in the two networks in at least one time step. 
The avalanche corresponding to a given knock-out is the set of perturbed genes (including 
the one which has been knocked-out), and its size is the number of perturbed genes. We have 
shown previously [Ml ES] that Saccharomyces cerevisiae shows a dynamical behaviour close 
to the “edge of chaos”, on the side of the ordered region and that for random topologies the 

^for the reasons discussed above 
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probability of having an avalanche involving m genes is given by the formula; 

( 3 ) 

where is a coefficienl|^ and A is equal to (1 — q){kout), where {kout) is the average of the 
number of outgoing connections per node. 

In [M] it was shown that RBN families with random topology, N = 6300 (the number 
of genes in Saccharomyces cerevisiae), kin = 2 and a particular set of Boolean functions (all 
the 13 two-input Boolean functions that remain after excluding the not biologically plausible 
XOR, EQUIVALENCE and FALSE functions). These choices allowed a very close ap¬ 
proximation of the experimental distributions. The measures we are proposing in this article, 
when applied to these networks, give respectively DA = 0.91 and SA = 0.94. These values 
of DA and SA are not very different, so if we were to use the former in Equation]^ we would 
obtain nonetheless a good approximation of the avalanche size distribution. In our previous 
papers on avalanche distribution [341 132) . we were interested in comparison with available 
experimental data. In this case no information is available about the features of the attractor 
state of the yeast, therefore it is correct to use average values such as DA and SA. However, 
in synthetic networks we can test the effect of the different attractor sensitivities on the re¬ 
sulting avalanches. By applying Equationj^to avalanches having different SAi, (i.e. different 
A’s) we obtain that the ratio between the frequency of avalanches of the same size m, starting 
from attractor a or attractor b, is ruled by the following equation: 

Note that these ratios are independent from the coefficients Bm. 

We checked this relationship by building 8000 RBN realizations, reaching one attractor for 
each realization (starting from random initial conditions), measuring its SAi value and per¬ 
forming a knock-out event. We also grouped the SAi values in bins, small but however 
sufficiently ample to have enough points in each bin. The comparison with the experimental 
results is good, as it can be seen in the three cases shown in Figure (corresponding to 
the fixed Xa values 0.94, 1.00 and 1.08). As shown in previous paragraphs, networks with 
a limited set of Boolean functions or evolved networks can have very different DA, SA and 
SAi'. in these cases the effect on the avalanche distributions are even more evident. A clear 
example is provided by A45 and AAq families: in these ensembles the theoretical and exper¬ 
imental DA are respectively equal to 0.75 and 0.74, but the SA measures are equal to 0.96 
and 0.67, these values originating from asymptotic states having respectively b = 0.08 and 
b = 0.67. These different asymptotic dynamical regimes give rise to very different avalanche 
distributions, as shown in figure Figure (note the different distribution tails). In this case 
the classical Derrida parameter DA would provide no information concerning the avalanche 
distribution. 


®See [23 for details on the calculation of these coefficients. 
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Figure 5: The comparison between the theoretical Rm ratios (computed from Equation Q 
and the experimental ones, the subset of results corresponding to m = 1 and maintaining 
fixed Xa respectively to 0.94, 1.00 and 1.08. The bars correspond to the standard deviation 
associated to each experimental Rm quotient: they are computed by propagating through 
the quotient the standard deviations associated to each experimental avalanche distribution, 
estimated by means of bootstrapping procedures. The experimental bins have a size of 0.01, 
each bin grouping from 110 till a maximum of 997 knock-out events. 
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Figure 6: Avalanche distribution in (x) and Ate (+) families. It is possible to observe 
the difference on the distribution tails. N = 6000, 10000 networks for each family. 
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7 Conclusion 


RBNs are widely used to simulate genetic regulatory systems, and it is well-known that they 
can show different dynamical regimes. In this paper we introduced different measures related 
to these dynamical behaviours, and we have shown that they provide a richer description of the 
dynamics than the usual ones, based on static sensitivity or on perturbations on random initial 
states. While some of the remarks presented here can be found also in other papers, which 
we have reviewed in this contribution, this work provides a more comprehensive analysis, and 
shows examples of a number of network families with somewhat unexpected properties. 

The main claims of this paper can be summarized as follows: 

• Like other dynamical systems, RBNs can show ordered, disordered and critical regimes. 
The same RBN can show several different asymptotic behaviours (attractors); 

• three kinds of measure are useful: (i) the classical Derrida measure (DA in the text), 
which characterizes each RBN structural parameters, (ii) the attractor sensitivity (S’vli 
in the text, relative to attractor i), which characterizes the dynamical regime of each 
attractor, and (Hi) the attractors sensitivity (SA in the text), which is related to the 
whole set of asymptotic dynamical behaviours of a RBN; 

• in classically critical RBNs DA and SA are very close, and they tend to coincide when 
considering ensemble averages; 

— single RBN realizations can show asymptotic behaviours with attractor sensitivity 
considerably far from its DA value; 

— and this tendency grows as the node connectivity kin increases; 

• finally, by knowing the asymptotic proportion of Os and Is it is possible to relate static 
and dynamic measures. 

We remark that the issues raised in this paper are not limited to RBNs, but they hold for 
a large class of discrete computational models and they can have strong consequences also for 
the biological processes that they describe. Indeed, these systems spend most of their time 
in their attractor cycles, so the most informative stability analysis should be performed in 
those regimes. The most common approach, based on perturbing random initial states, might 
fail in properly describing the asymptotic behaviour of these non-ergodic systems, except for 
some particular conditions. 
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